蒙特卡洛模拟¶
为什么需要MCMC¶
动机一¶
假如你需要对一维随机变量 $X$ 进行采样， $X$ 的样本空间是 $\{1,2,3\}$ ，且概率分别是 $\{1 / 2,1 / 4,1 / 4\}$ ，这很简单，只需写这样简单的程序：首先根据各离散取值的概率大小对 $[0,1]$ 区间进行等比例划分，如划分为 $[0,0.5],[0,5,0.75],[0.75,1]$ 这三个区间，再通过计算 机产生 $[0,1]$ 之间的伪随机数，根据伪随机数的落点即可完成一次采样。接下来，假如 $X$ 是连续 分布的呢，概率密度是 $f(X)$ ，那该如何进行采样呢? 聪明的你肯定会想到累积分布函数， $P(X<t)=\int_{-\infty}^t f(x) d x$ ，即在 $[0,1]$ 间随机生成一个数 $a$ ，然后求使得使 $P(x<t)=a$ 成立的 $t ， t$ 即可以视作从该分部中得到的一个采样结果。这里有两个前提: 一是概率密度函数可 积；第二个是累积分布函数有反函数。假如条件不成立怎么办呢? MCMC就登场了。
动机二¶
假如对于高维随机变量，比如 $\mathbb{R}^{50}$ ，若每一维取100个点，则总共要取 $10^{100}$ ，而已知宇宙的基本粒子大约有 $10^{87}$ 个，对连续的也同样如此。因此MCMC可以解决 "维数灾难" 问题。
均匀分布，Box-Muller 变换¶
均匀分布是很容易采样的，比如在计算机中生成 $[0,1]$ 之间的伪随机数序列，就可以看成是一种均匀分布。而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于 Uniform $(0,1)$ 的样本生成。例如正态分布可以通过著名的Box-Muller变换得到: Box-Muller变换 如果随机变量 $U_1, U_2$ 独立且 $U_1, U_2 \sim U n i$ form $[0,1]$ ， $$ \begin{aligned} &Z_0=\sqrt{-2 \ln U_1} \cos \left(2 \pi U_2\right) \\ &Z_1=\sqrt{-2 \ln U_1} \sin \left(2 \pi U_2\right) \end{aligned} $$ 则 $Z_0, Z_1$ 独立且服从标准正态分布。
其它几个著名的连续分布，包括指数分布、Gamma 分布、 $\mathrm{t}$ 分布、 $\mathrm{F}$ 分布、Beta 分布、Dirichlet 分布等等，都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。在python 的numpy， scikit-learn等类库中，都有生成这些常用分布样本的函数可以使用。 不过当 $p(x)$ 的形式很复杂，或者 $p(\mathbf{x})$ 是个高维的分布的时候，样本的生成就可能很困难了。譬如有如下的情况: 1) $p(x)=\frac{\tilde{p}(x)}{\int \tilde{p}(x) d x} ， \tilde{p}(x)$ 我们可以得到，但是下面的积分我们无法得到显示解（归一化系数 很难计算）； 2) $p(x, y)$ 是一个二维的分布函数，这个函数本身计算很困难，但是条件分布 $p(x \mid y), p(y \mid x)$ 的 计算相对简单; 如果 $p(\mathbf{x})$ 是高维的，这种情形就更加明显。
拒绝接受采样 (Acceptance-Rejection Sampling)¶
- 需求：根据已知分布的概率密度函数 $f(x)$ ，产生服从此分布的样本 $X$
- 准备工作:
- 需要一个辅助的 “建议分布proposal distribution" $G$ (已知其概率密度函数 $g(y)$ ) 来产生候选样本。由分布来产生候选样本，因此需要我们能够从 $G$ 抽样。即我们必须能够生成服从此概率分布的样本 $Y$ 。比如可以选择均匀分布、正态分布。
- 还需要另一个辅助的均匀分布 $U(0,1)$ 。
- 计算一个常数值 $c$ 。一一满足不等式 $c * g(x) \geqslant f(x)$ 的最小值 $c$ (当然，我们非常希望 $c$ 接近于1)
- 开始样本生成:
- 从建议分布 $G$ 抽样，得到样本 $Y$ 。
- 从分布 $U(0,1)$ 抽样，得到样本 $U$ 。
- 如果 $U \leqslant \frac{f(Y)}{c * g(Y)}$ ，则令 $X=Y \quad($ 接受 $Y)$ ，否则继续执行步骤1 (拒绝)。
接受拒绝采样的直观解释¶
假设使用 $U(0,1)$ 来作为 “proposal distribution" $G$ ，这样 $g(x)=1 \forall x \in[0,1]$ 。如下图所示，我们每次生成的两个样本 $Y$ 与 $U$ ，对应下图中矩形内的一点 $P(Y, U * c * g(Y))$ 。接受条件 $U \leqslant f(Y) c * g(Y)$ ，即 $U * c * g(Y) \leqslant f(Y)$ 的几何意义是点 $P$ 在 $f(x)$ 下方，不 接受 $Y$ 的几何意义是点 $P$ 在 $f(x)$ 的上方。在 $f(x)$ 下方的点(o形状)满足接受条件，上方的点 (+形状)不满足接受条件。
python实现¶
假如我们的目标概率密度函数是 $f(x)=\frac{(x-0.4)^4}{\int_0^1(x-0.4)^4 d x}$ ，对此分布生成样本。
- 准备工作:
- 建议分布函数 $G$ 选择 $U(0,1)$ ，概率密度函数为 $g(x)=1 \forall x \in[0,1]$ 。（这里我们简单化 了，一般要使得c最小的g函数)；
- 辅助的均匀分布 $U(0,1)$ ；
- 计算常数值， $c=\max \left\{0.4^4 / \mathrm{sum}, 0.6^4 / \mathrm{sum}\right\}=0.6^4 / \mathrm{sum}$ ，(其中 $s u m=\int_0^1(x-0.4)^4 d x ，$ 后面可以化简掉)
import random
import math
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
%matplotlib inline
sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = (12, 8)
def AceeptReject(split_val):
global c
global power
while True:
x = random.uniform(0, 1)
y = random.uniform(0, 1)
if y*c <= math.pow(x - split_val, power):
return x
power = 4
t = 0.4
sum_ = (math.pow(1-t, power + 1) - math.pow(-t, power + 1)) / (power + 1) #求积分
x = np.linspace(0, 1, 100)
#常数值c
c = 0.6**4/sum_
cc = [c for xi in x]
plt.plot(x, cc, '--',label='c*f(x)')
#目标概率密度函数的值f(x)
y = [math.pow(xi - t, power)/sum_ for xi in x]
plt.plot(x, y,label='f(x)')
#采样10000个点
samples = []
for i in range(10000):
samples.append(AceeptReject(t))
plt.hist(samples, bins=50, normed=True,label='sampling')
plt.legend()
plt.show()
蒙特卡罗方法小结¶
使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和得到。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如之前的第一个问题有时可以解决，但又会产生另一个问题:
- 对于一些二维分布 $p(x, y)$ ，我们只能得到条件分布 $p(x \mid y), p(y \mid x)$ ，却不能得到二维分布的 一般形式；
- 对于复杂非常见分布 $p\left(x_1, x_2, \ldots, x_n\right)$ ，我们很难找到合适的 $q(x), c$ 。 要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，还的需马尔科夫链的帮忙。
马尔科夫采样¶
马尔科夫链的收敛性质¶
如果一个非周期的马尔科夫有状态转移矩阵 $P$ ，并且他的任何两个状态是连通的，那么 $\lim _{n \rightarrow \infty} P_{i j}^n$ 与 $i$ 无关，我们有:
- $\lim _{n \rightarrow \infty} P_{i j}^n=\pi(j)$
- $\lim _{n \rightarrow \infty} P^n=\left[\begin{array}{ccccc}\pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \ldots \\ \vdots & \vdots & \vdots & \ddots & \vdots\end{array}\right]$
- $\pi$ 是方程 $\pi P=\pi$ 唯一非负解，其中 $\pi=[\pi(1), \pi(2), \cdots, \pi(j), \cdots] \sum_{i=1}^{\infty} \pi(i)=1$ ， $\pi$ 称为马氏链的平稳分布。
基于马尔科夫链采样¶
如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采用出这个平稳分布的样本集。 首先，基于初始任意简单概率分布比如高斯分布 $\pi_0(x)$ 采样得到状态值 $x_0$ ，基于条件概率分布 $P\left(x \mid x_0\right)$ 采样状态值 $x_1$ ，一直进行下去，当状态转移进行到一定的次数时，比如到 $n$ 次时，我 们认为此时的采样集 $\left(x_n, x_{n+1}, \ldots\right)$ 即是符合我们的平稳分布的对应样本集，可以用来做蒙特 卡罗模拟求和了。 算法如下：
- 输入马尔科夫链状态转移矩阵 $P$ ，设定状态转移次数阈值 $n_1$ ，需要的样本个数 $n_2$ ；
- 从任意简单概率分布采样得到初始状态值 $x_0$ ；
- for $t=0$ to $n_1+n_2-1$ : 从条件概率分布 $P\left(x \mid x_t\right)$ 中采样得到样本 $x_{t+1}$ ；
- 样本集 $\left(x_n, x_{n+1}, \ldots\right)$ 即是符合我们的平稳分布的对应样本集。
马尔科夫链采样小结¶
如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵，那么我们 就可以用马尔科夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。 但是一个重要的问题是，随意给定一个平稳分布 $\pi$ ，如何得到它所对应的马尔科夫链状态转移矩阵 P呢? MCMC是时候出现了。
马尔可夫链蒙特卡罗算法¶
马尔科夫链的细致平稳条件(Detailed Balance Condition)¶
在解决从平稳分布 $\pi$ ，找到对应的马尔科夫链状态转移矩阵 $P$ 之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下: 如果非周期马尔科夫链状态转移矩阵 $P$ 和概率分布 $\pi(x)$ 对所有的 $i, j$ 满足: $$ \pi(i) P(i, j)=\pi(j) P(j, i), \text { for all } i, j $$ 则称概率分布 $\pi(x)$ 是状态转移矩阵 $P$ 的平稳分布(Stationary Distribution)。 上述只是个充分条件，证明如下: 因为 $\sum_i \pi(i) P(i, j)=\sum_i \pi(j) P(j, i)=\pi(j) \sum_i P(j, i)=\pi(j)$ ，所以 $\pi P=\pi$ 。 当分布是二维时，此条件是充要的，但3维以上时，就不是了。
从细致平稳条件可以得到，只要我们找到了可以使概率分布 $\pi(x)$ 满足细致平稳分布的矩阵 $P$ 即 可。这给了我们寻找从平稳分布 $\pi$ ，找到对应的马尔科夫链状态转移矩阵 $P$ 的新思路。 不过不幸的是，我们很难找到合适的矩阵 $P$ 满足细致平稳条件。
MCMC采样¶
由于一般情况下，目标平稳分布 $\pi(x)$ 和某一个马尔科夫链状态转移矩阵 $Q$ 不满足细致平稳条 件，即： $$ \pi(i) Q(i, j) \neq \pi(j) Q(j, i) $$ 以下记号表达同一个意思: $Q(i, j) \Leftrightarrow Q(j \mid i)) \Leftrightarrow Q(i \rightarrow j)$ 我们引入一个 $\alpha(i, j)$,使上式可以取等号。 $$ \pi(i) Q(i, j) \alpha(i, j)=\pi(j) Q(j, i) \alpha(j, i) $$ 怎样才能成立呢，其实很简单，按照对称性： $$ \alpha(i, j)=\pi(j) Q(j, i) ; \pi(i) Q(i, j)=\alpha(j, i) $$ 然后我们就可以得到了分布 $\pi(x)$ 对让马尔科夫链状态转移矩阵 $P(i, j)=Q(i, j) \alpha(i, j)$ 其中 $\alpha(i, j)$ 一般称之为接受率，取值在 $[0,1]$ 之间，可以理解为一个概率值。这很像接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布， 这里是以一个常见 的马尔科夫链状态转移矩阵 $Q$ 通过一定的接受-拒绝概率得到目标转移矩阵 $P$ ，两者的解决问题思路是类似的。
MCMC采样算法如下：
- 输入任意给定的马尔科夫链状态转移矩阵 $Q$ ，目标平稳分布 $\pi(x)$ ，设定状态转移次数阈值 $n_1$ ，需要的样本数 $n_2$ ；
- 从任意简单概率分布得到初始状态值 $x_0$ ；
- for $t=0$ in $n_1+n_2-1$ a. 从条件概率分布 $Q\left(x \mid x_t\right)$ 得到样本值 $x_*$ b. 从均匀分布中采样 $U \sim[0,1]$ c. 如果 $u<\alpha\left(x_t, x_*\right)=\pi\left(x_*\right) Q\left(x_*, x_t\right)$ ，则接受 $x_t \rightarrow x_*$ ，即 $x_{t+1}=x_*$ d. 否则不接受转移， $t=\max \{t-1,0\}$ 上述过程中 $p(x), q(x \mid y)$ 说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然 是有效，于是就得到更一般的连续概率分布 $p(x)$ 的采样算法，而 $q(x \mid y)$ 就是任意一个连续二元 概率分布对应的条件分布。 但是这个采样算法还是比较难在实际中应用，因为在第三步中，由于 $\alpha\left(x_t, x_*\right)$ 可能非常的小， 比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个 $n_1$ 要非常非常的大，这让人难以接受，怎么办呢? 这时就 轮到我们的M-H采样出场了。
M-H采样¶
M-H采样算法¶
M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样。 M-H采样解决了我们上一节MCMC采样接受率过低的问题。 我们可以对 $\pi(i) Q(i, j) \alpha(i, j)=\pi(j) Q(j, i) \alpha(j, i)$ 两边进行扩大，此时细致平稳条件也是 满足的，我们将等式扩大 $C$ 倍，使 $C \alpha(j, i)=1$ (精确来说是使得两边最大的扩大为 1 )，这 样我们就提高了采样中的跳转接受率，所以我们可以取 $\alpha(i, j)=\min \left\{\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}, 1\right\}$
M-H采样算法如下:
- 输入任意给定的马尔科夫链状态转移矩阵 $Q$ ，目标平稳分布 $\pi(x)$ ，设定状态转移次数阈值 $n_1$ ，需要的样本数 $n_2$ ；
- 从任意简单概率分布得到初始状态值 $x_0$ ；
- for $t=0$ in $n_1+n_2-1$ a. 从条件概率分布 $Q\left(x \mid x_t\right)$ 得到样本值 $x_*$ b. 从均匀分布中采样 $U \sim[0,1]$ c. 如果 $u<\alpha\left(x_t, x_*\right)=\min \left\{\frac{\pi(*) Q\left(x_*, x_t\right)}{\pi(t) Q\left(x_t, x_*\right)}, 1\right\}$ ，则接受 $x_t \rightarrow x_*$ ， 即 $x_{t+1}=x_*$ d. 否则不接受转移， $t=\max \{t-1,0\}$
很多时候我们马尔科夫转移状态矩阵是对称的，因此接受率可以写成: $\alpha(i, j)=\min \left\{\frac{\pi(j)}{\pi(i)}, 1\right\}$ 对于分布 $p(x)$ ，我们构造转移矩阵 $P$ 使其满足细致平稳条件: $$ \pi(x) P(x \rightarrow y)=\pi(y) P(y \rightarrow x) $$ 此处 $x$ 并不要求是一维的，对于高维空间的 $p(\mathbf{x})$ ，如果满足细致平稳条件 $$ \pi(\mathbf{x}) P(\mathbf{x} \rightarrow \mathbf{y})=\pi(\mathbf{y}) P(\mathbf{y} \rightarrow \mathbf{x}) $$ 那么以上的 Metropolis-Hastings 算法一样有效。
M-H采样python实现¶
假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 $Q(i, j)$ 的条件转移概率是以 $i$ 为均值，方差 1 的正态分布在位置 $j$ 的值。
from scipy.stats import norm
def norm_dist_prob(theta):
y = norm.pdf(theta, loc=3, scale=2)
return y
T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
t = t + 1
pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None) #状态转移进行随机抽样
alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1]))) #alpha值
u = random.uniform(0, 1)
if u < alpha:
pi[t] = pi_star[0]
else:
pi[t] = pi[t - 1]
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2),label='Target Distribution')
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')
plt.legend()
plt.show()
M-H采样小结¶
在大数据时代， M-H采样面临着两大难题: 1) 我们的数据特征非常的多，M-H采样由于接受率计算式 $\min \left\{\frac{\pi(j) Q(j, i)}{\pi(i) Q(i, j)}, 1\right\}$ 的存在，在高维 时需要的计算时间非常的可观，算法效率很低。同时 $\alpha(i, j)$ 一般小于 1 ，有时候辛苦计算出来却 被拒绝了。能不能做到不拒绝转移呢? 2) 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出 各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采 样呢?
Gibbs采样¶
二维Gibbs采样¶
利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。算法如下:
- 输入平稳分布 $\pi\left(x_1, x_2\right)$ ，设定状态转移次数阈值 $n_1$ ，需要的样本个数 $n_2$ ；
- 随机初始化初始状态值 $x_1^{(1)}$ 和 $x_2^{(1)}$ ；
- for $t=0$ to $n_1+n_2-1$ : -从条件概率分布 $p\left(x_2 \mid x_1^t\right)$ 中采样得到样本 $x_2^{t+1}$
- 从条件概率分布 $p\left(x_1 \mid x_2^{t+1}\right)$ 中采样得到样本 $x_1^{t+1}$ 样本集: $\left\{\left(x_1^{\left(n_1\right)}, x_2^{\left(n_1\right)}\right),\left(x_1^{\left(n_1+1\right)}, x_2^{\left(n_1+1\right)}\right), \ldots,\left(x_1^{\left(n_1+n_2-1\right)}, x_2^{\left(n_1+n_2-1\right)}\right)\right\}$ 即为我们需要 的平稳分布对应的样本集。 整个采样过程中，马氏链的转移只是轮换的沿着坐标轴，采样的过程为: $$ \left(x_1^{(1)}, x_2^{(1)}\right) \rightarrow\left(x_1^{(1)}, x_2^{(2)}\right) \rightarrow\left(x_1^{(2)}, x_2^{(2)}\right) \rightarrow \cdots \rightarrow\left(x_1^{\left(n_1+n_2-1\right)}, x_2^{\left(n_1+n_2-1\right)}\right) $$
马氏链收玫后，最终得到的样本就是 $p(x, y)$ 的样本，而收玫之前的阶段称为 burn-in period。 额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其 实是不强制要求的。最一般的情形可以是，在 $t$ 时刻，可以在 $x$ 轴和 $y$ 轴之间随机的选一个坐标 轴，然后按条件概率做转移，马氏链也是一样收玫的。轮换两个坐标轴只是一种方便的形式。 $6.3$ 多维Gibbs采样 上面的这个算法推广到多维的时候也是成立的。比如一个 $n$ 维的概率分布 $\pi\left(x_1, x_2, \ldots, x_n\right)$ ， 可以通过在 $n$ 个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴 $x_i$ 上的转 移，马尔科夫链的状态转移概率为 $P\left(x_i \mid x_1, x_2, \ldots, x_{i-1}, x_{i+1}, \ldots, x_n\right)$ ，即固定 $n-1$ 个 坐标轴，在某一个坐标轴上移动，即每次转移只允许某一个维度上变化，可以解决维数灾难问题。 算法如下:
- 输入平稳分布 $\pi\left(x_1, x_2, \ldots, x_n\right)$ ，设定状态转移次数阈值 $n_1$ ，需要的样本个数 $n_2$
- 随机初始化初始状态值 $\left.x_1^{(1)}, x_2^{(1)}, \ldots, x_n^{(1)}\right)$
- for $t=0$ to $n_1+n_2-1$ :
- 从条件概率分布 $p\left(x_1 \mid x_2^t, x_3^t, \ldots, x_n^t\right)$ 中采样得到样本 $x_1^{t+1}$
- 从条件概率分布 $p\left(x_2 \mid x_1^{t+1}, x_3^t, \ldots, x_n^t\right)$ 中采样得到样本 $x_2^{t+1}$
- ....
- 从条件概率分布 $p\left(x_j \mid x_1^{t+1}, x_2^{t+1}, x_{j-1}^{t+1}, x_{j+1}^t \ldots, x_n^t\right)$ 中采样得到样本 $x_j^{t+1}$
- ....
- 从条件概率分布 $p\left(x_n \mid x_1^{t+1}, x_2^{t+1}, \ldots, x_{n-1}^{t+1}\right)$ 中采样得到样本 $x_n^{t+1}$
样本集: $\left\{\left(x_1^{\left(n_1\right)}, x_2^{\left(n_1\right)}, \ldots, x_n^{\left(n_1\right)}\right), \ldots,\left(x_1^{\left(n_1+n_2-1\right)}, x_2^{\left(n_1+n_2-1\right)}, \ldots, x_n^{\left(n_1+n_2-1\right)}\right)\right\}$ 即为 我们需要的平稳分布对应的样本集。当然这些样本并不独立，但是我们此处要求的是采样得到的样 本符合给定的概率分布，并不要求独立。
整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定 $n-1$ 个特 征，对某一个特征求极值。而Gibbs采样是固定 $n-1$ 个特征在某一个特征采样。
同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的 Gibbs采样的实现都是基于坐标轴轮换的。
二维Gibbs采样实例python实现¶
假设我们要采样的是一个二维正态分布 $N(\mu, \Sigma)$ ，其中： $\mu=\left(\mu_1, \mu_2\right)=(5,-1)$ ， $\Sigma=\left(\begin{array}{cc}\sigma_1^2 & \rho \sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma_2^2\end{array}\right)=\left(\begin{array}{ll}1 & 1 \\ 1 & 4\end{array}\right)$ 而采样过程中的需要的状态转移条件分布为: $$ P\left(x_1 \mid x_2\right)=N\left(\mu_1+\rho \sigma_1 / \sigma_2\left(x_2-\mu_2\right),\left(1-\rho^2\right) \sigma_1^2\right) $$ $$ P\left(x_2 \mid x_1\right)=N\left(\mu_2+\rho \sigma_2 / \sigma_1\left(x_1-\mu_1\right),\left(1-\rho^2\right) \sigma_2^2\right) $$
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])
def p_ygivenx(x, m1, m2, s1, s2):
return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))
def p_xgiveny(y, m1, m2, s1, s2):
return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))
N = 5000
K = 20
x_res = []
y_res = []
z_res = []
m1 = 5
m2 = -1
s1 = 1
s2 = 2
rho = 0.5
y = m2
for i in range(N):
for j in range(K):
x = p_xgiveny(y, m1, m2, s1, s2) #y给定得到x的采样
y = p_ygivenx(x, m1, m2, s1, s2) #x给定得到y的采样
z = samplesource.pdf([x,y])
x_res.append(x)
y_res.append(y)
z_res.append(z)
num_bins = 50
plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5,label='x')
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5,label='y')
plt.title('Histogram')
plt.legend()
plt.show()
Gibbs采样小结¶
由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当 然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一 维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。 有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。
资料来源
